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ABSTRACT 





The ability to measure and compensate for power line harmonics has become a 


growing area of interest because of today’s commonly used electronic equipment. Since 


the number and relative magnitudes of the harmonics on the power line are a function 


of the load, estimation of an equivalent load can be accomplished. Because of variation 


in the load, the need for an adaptive algorithm is imperative. Thus far, few algorithms 


for determining harmonic contents have not dealt with the problem associated with the 


limited power of the line conditioner. 


This thesis deals with a previously known harmonic compensating algorithm and 


introduces a new algorithm which both compensates for harmonic noise and estimates the 


load as a transformation matrix depending on the associated transfer function of the 


active line conditioner in use. 
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I. INTRODUCTION 


A. BACKGROUND 

Many of today’s electronic devices i.e., computers, fluorescent lights and 
microwave ovens, effect power distribution due to their nonlinear consumption of power. 
The result is irregular current and voltage wave forms on the power line. Active line 
conditioners provide a way of eliminating the accompanying noise on a power line by 
independently adjusting the active and nonactive components, thereby maintaining a 


constant sinusoidal bus voltage. 


line 
conditioner 


nonactive 


Figure 1.1 Bus with Linear and Nonlinear Loads and Active Line Conditioner 


The amount of voltage distortion is a function of the nonlinear load distribution, 
their current spectra, topology and frequency dependance within the network. For a 


nonsinusoidal voltage [Ref 1] 


V, = v2V,sinot + /2)° V,sin(hot + «,) 
hel 





with linear and non linear loads drawing a total current. 





=i thy + by (1.2) 
where 
i, = ¥21,cos®, sinot (1.3) 
is the active component of the fundamental current, 
i, = ¥21,sin®,cosut (1.4) 
is the reactive component of the fundamental current and 
ig = v2) I,sin(hot + «, + 9,) (1.5) 


hel 


«€ 
is the harmonic current. The apparent power can be described as follows. 


S=ViI, = G *+Pwoag +O (1.6) 
hel hel 


Equation 1.6 can be expressed in phasoz form, 


$= |, + PY + + Oy oe) 
where 
P, = V,1,cos6, (1.8) 


represents the Fundamental Power Frequency Active Power. The Harmonic Active 


Power is, 











P, = > V,1,c0s8, (1.9) 


and 
Q, = V,J,sin6, (1.10) 


is the Power Frequency Reactive Power. The Harmonic Reactive Power is 


Q,y = > V,1,sin6, (1.11) 


hel 


The components of Q,, are generated by specific harmonic voltages and harmonic 
currents (in no particular order) [Ref 2]. Because of the vector properties of the reactive 
power components, control or cancellation is possible using vectors of identical frequency 
and magnitude with opposite phase. Therefore, the Harmonic Reactive Power can be 
eliminated by introducing or drawing current from the power line which is 180° out of 
phase with each respective harmonic. With this in mind, the conditioner current is as 


follows, 


ig ig tie the (1.12) 


where ic,, and ic,, are the conditioner Fundamental Active and Reactive currents 


respectively, and 


iy = ee ee +O, + Y,) (1.13) 


is the conditioner harmonic current. 


B. ACTIVE LINE CONDITIONERS 


Active line conditioners serve a dual purpose. First, they adjust one or more loads 


thereby changing the active power. Secondly, they are capable of controlling the 


amplitude and phase characteristics of the nonreactive currents i,,, and i,,, which affect 
the value of Q, and/or Q,, [Ref 2,4]. By using a solid state switching network at a 
frequency much greater then that of the fundamental, the line current is modulated in 
order to maintain the boundaries of a desired template z, shown in Figure 1.2. Fora 
narrow boundary of error 6, the conditioner current i, is unaffected by the fluctuations 
(the zig-zaging) within the boundary area. By adjusting the template waveform through 
a feedback circuit, the conditioner current spectrum can be altered to produce an overall 
current I, that is as close to sinusoidal as possible. To do this, the load current i, is 
monitored to determine the harmonics present. Then by injecting a current i, from the 
line conditioner which is 180° out of phase with that of the load, the unwanted harmonics 


can be cancelled out. 


1. Equivalent Circuit Modeling 
For medium and low voltage systems, the best practical means of adjusting 
the conditioner current i, is by minimization of the voltage distortions at the conditioner 


location on the bus [Ref. 5]. The voltage at the conditioner node can be represented by 


N 
Ven = 2 Zarb (1.14) 
pas 


where 


I,, = Bus x harmonic current phasor of order h 











template z 





0 Tv 


Figure 1.2 Bus Voltage, Line Conditioner Template and Current 


Zan, = Harmonic Complex Impedances of the node 
N = Number of independent nodes 


A Norton equivalent circuit for the bus has the equivalent harmonic current 


N 
1, = Yeas » Zorvlin (1.15) 
xel xeB 


where Yp—, = 1/Zgg, is the self-admittance of the node for the harmonic h. Figure 1.3 


shows the Norton equivalent circuit with the associated harmonic currents and load Z,,. 


v Deo Deeg , 


Figure 1.3 System Approximation using Norton Equivalent Circuit 


Where Zy, is given in Equation 1.16, 





Zs, = Roy + fhols, = Zyp,1Zp, (1.16) 


and Z,, is the load impedance of the harmonic order h. 
From the equivalent load Z,, the voltage due to the offending harmonics V,, 


can be defined as 
Vin = Len Zsp (1.17) 


For a linear resistance and inductor Rg, and L,,, a conditioner current i. equal 
to the negative of I,, would eliminate the harmonic voltage V,,,. It is important to note 
that active line conditioners are inherently limited in their maximum current output, 
therefore, negating the entire value of I,, may not be possible. Although limited, any 
reduction of the harmonic noise, especially of lower order, significantly improves the 


recognition of the fundamental. 


II. SURVEY OF PERVIOUS WORK 


A. ADAPTIVE ESTIMATION OF HARMONIC VOLTAGE 
The best fitting sinusoidal wave to a nonsinusoidal periodic wave is the fundamental 


{Ref. 1]. The error associated with such a system can be written. 
€ = Vp - VY = Vp + Veo (2.1) 


Since the signals v, and v, are periodic, the error, ¢, is Statistically stationary. 
Therefore, the expected value of the square of the error ¢€ results in a quadratic function 
which has a guaranteed minimum for real physical signals (Ref. 6]. Then by minimizing 
the mean square of the error (MSE), the signal should be nearly identical to that of the 
fundamental. By representing the error voltage due to the harmonics as the sum of 
weighted sines and cosines, an error surface for each weight can be defined thereby 
making the calculation of a minimum possible. Figure 2.1 shows the block diagram for 
such an adaptive system. 

Xg2, Xs3 “- Xsn aNd X~, Xc3 ~* Xcp Tepresent discrete versions of the harmonics 
associated with the power line, and C,, C;, -- C, and S,, S;, -- S, are the weights of the 


associated harmonic. The correction voltage to the conditioner is defined by 


y, = ¥ ([S,sin(hw@kT/N) + C,cos(hokT/N)) (2.2) 


hal 


where k = time index 
T = 1/f = 2x/w = period or one cycle 





Figure 2.1 Block Diagram of Adaptive System 


N = number of samples per cycle 


The Active Line Conditioner current is given by 


icy = ae > [S,sin(hwt) + C,cos(hwt)] 


where D = The gain of the D/A converter 
K = Converter constant of the Conditioner in (A/V) 


The discrete error e, as a function of voltage becomes 


ey = Vox * Vi 


where Vg, and v,, are represented by 


Yee = KDY Z,{5,8in(hokTIN) + C,cos(hwkT]N)) 


(2.3) 


(2.4) 


(2.5) 





Vin = 2 Z,1,,[cosB,sin(hokT/N)+ sinB,cos(hwkT/N)) (2.6) 


Simply setting the sine and cosine weights equal to 


S, = -v2U,,cosB,)/KD (2.7) 
C, = -¥2(,,c0sB,)/KD 


requires the error to be zero. Since the actual load impedance Z,, is not known and 
changes with time, an estimation of the sine and cosine weights is performed by the MSE 


processor using the following linear prediction. 


Spat iz 5, * (-V,) n/h 


(2.8) 
Cyr = 5 + (“Ved ala 
where 
V, = de/os 
2.9 
Vo = de/dCc os 


are the error gradients of the sine and cosine weights respectively, and yu is a constant 
called the acceleration factor which is directly related to both the rate of convergence, 
and the magnitude of any over-shoot in reaching the minimum. The / term is used to 
scale the acceleration factor by an amount proportional to the harmonic being evaluated. 
This allows for faster convergence of lower order harmonics, the largest error, without 
driving the higher order unstable. Figure 2.2 shows a quadratic error surface as a 


function of a single weight. 





MEAN SQUARE ERROR ¢é 


Cyr4 Ce +3 
WEIGHT 


Figure 2.2 Single Weight Error Surface Example 


Starting on the left side where C,,, > C, and &4, < & the gradient is negative 
and the value of ¢ converges toward the minimum. If the value of the weighting factor 
produces a higher error then the previous &,, > &, 3» an Over shoot occurs, but the 
gradient remains negative thereby predicting a smaller weighting value then the current 
one and € again converges toward the minimum. The same result would be obtained for 


an initial weight greater then that need to minunize &. 


B. LINEAR LOAD SIMULATION 


A program for testing the validity of the adaptive algorithm was written for 


MATLAB using the parameters in Figure 2.3 and a noise component equating to ten 


percent Total Harmonic Distortion (THD) assuming only odd harmonics up to the 
twenty-first [Ref. 1]. The fundamental frequency is assumed to be a standard 60 Hz 


outlet and the line conditioner to have no power limitation. 








Figure 2.3 Linear Load Equivalent Circuit Model 


Figure 2.4 shows the MSE for the first three offending harmonics and their 
convergence to zero, all well under one second using an acceleration constant of 3e-7. 
The THD in Figure 2.5 also follows a similar pattern since it is most affected by the 
lower order harmonics and becomes one-one hundredth of its original value after just one 
half second, or thirty iterations. Figure 2.6 shows the weighting coefficients of the sine 
and cosine for the first three offending harmonics. With the help of a simple 
trigonometric identity, these values can easily be shown to correspond to the magnitude 
and phase of the original noise components. Also note that the final weighting values 
where reached asymptotically without any over-shoot indicating the choice of the 
acceleration coefficient is optimal with respect to requiring minimum power from the 


active line conditioner. 
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Figure 2.4 Mean Square Error for Linear Load 
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Total Harmonic Distortion 
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Figure 2.5 Total Harmonic Distortion Linear Load 
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Figure 2.6 Sine & Cosine Weighting Coefficients for Linear Load 
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C. NONLINEAR LOAD SIMULATION 


The inductor in Figure 2.3 was substituted for the one shown in Figure 2.7 to 





produce a nonlinear response in the load. The nonlinearity was chosen to provide 


approximately ten percent deviation from that of the linear case. 





Figure 2.7 Nonlinear Load Equivalent Circuit Model 


The nonlinear load results were very similar to those from the linear with small 
deviations in the THD and identical results for the weighting coefficients. Some of these 
values are summarized below in Table 2.1. 


TABLE 2.1 TOTAL HARMONIC DISTORTION FOR LINEAR AND NONLINEAR 
LOADS 










| 
{ 
! 


[Linear | 9.327 | 8.89 | 6.381 | 3.300 | 1.067 | 123 | .oo1 | 
[Nonlinear | 9.245 | 3.874 | 6.410 | 3.427 | 1.120 | 120 | 0012 | 
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Mean Square Error 
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Figure 2.8 Mean Square Error Nonlinear Load 
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Figure 2.9 Total Harmonic Distortion Nonlinear Load 
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Figure 2.10 Sine & Cosine Weighting Coefficients for Nonlinear Load 
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TI. OPTIMAL ESTIMATION WITH SYSTEM IDENTIFICATION 


A. MODELING THE NETWORK 

The model used for optimal estimation is very similar to that found in Figure 2.1 
with the exception that the impedance, although unknown, will be estimated along with 
minimizing the harmonic error. The block diagram model is shown in Figure 3.1 with 
the impedance of the present load represented as a square matrix H. It should be noted 
that this model can be used for different system transfer function input/output 


relationships other than current to voltage. 


w(n) 


u(n) H yn) (+) e(n) 


Figure 3.1 Optimal Estimation Impedance Model 


Since the harmonic noise w(t) is assumed to be harmonics of the fundamental it can 


be represented as a sum of sinusoids shown in Equation 3.1 for a continuous signal. 
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Ny 


w(t) = }> A,cos(2xxf,) + B,sin(2xxf,0) (3.1) 


x=2 


For a discrete sampled system Nyquist criteria must be maintained, therefore the 
sampling frequency, /,, must be an integer value of the fundamental which is greater then 
twice the highest harmonic frequency to be eliminated. The continuous frequencies are 


converted to discrete under the following conditions. 


ant = ant where 1<1<N,< M (3.2) 
f M 2 


RY 


From this the disturbance w(t) can be represented discretely in matrix form as 


w(n) = W'x(n) (3.3) 


w' = [0, Q, A,, B,, A;, B;, Tore Ay,» By] 


x(n) = [cos2x Ln, sin2x Ln, cos2n 2m, sin 2x 2-2, face 5 
M M M M 


N, N 
cos2n—n, sin2x—tn]™ 
M M 


Because the harmonic noise does not change instantaneously with changes to the 
load, it is reasonable to assume that it is periodic. By dividing the time scale into 


periodic intervals of length N, which is a multiple of M, then for all n x(n) = x(” + M) 


= x(n + N). 











HHH HHH 


N (k + 1)N (k + 2)N 








Figure 3.2 Discrete Time Period 






By defining the control input current to be 






un) = usx(n) KN sus KN+N-1 (3.5). 


where u, is a weight vector to be determined. The error can now be written discretely 







in terms of x(n) as 







e(n) = W'x(n) + u, H' x(n) 






Note that for a linear impedance H becomes a diagonal matrix, otherwise it is not. 


Since w(n) is the signal which is to be eliminated, by using its frequency 







information in the control input, u(m), it will provide some of the needed information to 






negate it. The remaining control information will come from a recursive estimation of 


its Fourier coefficients which is the main basis of this algorithm. 










B. CONTROL 






Since the value of H can be estimated through the use of system identification 


techniques, it can be assumed to be known for the purpose of determining the control 







necessary to eliminate the harmonic noise. From Equation 3.6 it is easy to see that W" 





and u,"H" are scalars which can be combined to represent the error weights for each 


sample point over a single interval, N, of the fundamental. The error associated with 


a single sample n is shown in Equation 3.7. 


e(n) = e, x(n) 
Now each term in Equation 3.6 can be defined in terms of x(n). 
Due to the periodicity of the harmonic frequencies, the frequency components can 
be eliminated by subtracting the error of the k” interval from the (k-1)" interval resulting 


in a difference equation of just weighted vectors shown below. 


e, x(n) = W'x(n) + u,H™x(n) 


e,-,X(2-N) = W'x(n-N) + u,.,H™x(n-N) 
e;- &-. = AT; - u-;) 


After some manipulation, Equation 3.8 can be written as 


Qle, - &_,) = Uy ~ %., 
where Q = H' 
Now that the control is in terms of the error and an admittance matrix Q, using a 


linear predictor similar to that in Equation 2.8 can be used. 


where a is a scalar defined on the interval -1< a <1. 


The error, e,, can be driven to zero for Q not equal to the null set. 








C. ESTIMATION 


Because the load to the bus changes over time some method of updating the value 
of H to those changes must exist in order for the controller to effectively reduce the 
harmonic noise present. Estimation of the admittance matrix Q can easily be 
incorporated with the control through system identification methods using a recursive 


least squares (RLS) algorithm. By choosing estimates of the control and error to be 


iT} 
Pw 
i 


=u, ~ UY) 


: (3.11) 
ee, — &, 
Equation 3.9 reduces to a matrix form of the RLS equation. 
a, = &Q (3.12) 


Although Equation 3.12 is in matrix form, it is important to remember that the 
estimate of each row of Q is a unique difference equation of the associated control and 
error coefficients of their respective frequency. In other words, even though the 
frequency vector x(m) is not found in Equation 3.12 the relationship between the third 
harmonic error and control coefficients remains linear. It is well known that the output 
of a linear system differs from the input only in magnitude and phase, therefore the 


system output y(t) can be written in terms of u(t) as follows: 


pe) 


y,(t) = |H,|A,cos(2nkfr+a,) + |H,|B,sin(2nkft+a,) = 


C.cos(2nkf,t) + Dysin(2wkf,t) Oe 


where C, and D, reflect the magnitude and phase changes of the system on the input 
coefficients A, and B,. With the help of a trigonometric identity, Equation 3.13 can be 


written in a more convenient matrix form [Ref 7}. 


“1H, 


Since H,, cosa, and sina, are scalars Equation 3.14 can be arrange in a recursive form 


cosa, -sina, 
sina, cosa, 











i. (3.14) 
B, 





similar to that of Equation 3.12 as follows 





ct - au Mi (3.15) 
Dy} [Be Ax) [2x 
where Y, and Z, represent the estimate for the product of the magnitude and phase 
change for the cosine and sine terms respectively of a given harmonic. Equation 3.15 
represents the estimates of the coefficients for just a single harmonic of the system and 
can be thought of as a building block of the matrix RLS Equation 3.12. 

Now all the needed information is available for implementation of a Kalman filter 
based RLS estimation [Ref 8, 9]. Figure 3.3 shows a block diagram representation of 


the system model with the estimation algorithm. 


1. Considerations in Applying FFT to Harmonic Analysis 
Since this algorithm emphasizes the use of Fourier coefficients, some aspects 


of using the Fast Fourier Transform (FFT) will be addressed. The FFT algorithm has 
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ESTIMATOR 





Figure 3.3 System Diagram with Estimation and Control 


useful applications in power system networks, but can produce erroneous information if 
not applied correctly. Certain assumptions about the FFT must be understood to avoid 
false representation of the associated signal [Ref 7]. 

© The signal is stationary (constant magnitude). 

@ Each frequency in the signal is an integer multiple of the fundamental. 


@ The sampling frequency is equal to the number of samples times the fundamental 
frequency used in the algorithm. 


@ The sampling frequency meets Nyquist criteria. 


D. SIMULATION RESULTS 

In testing the optimal estimation algorithm the same harmonic noise components 
from the MSE in Chapter II were used. Since the impedance matrix of the system is 
estimated using RLS any linear transformation matrix for H can be used. For the 


purpose of this research a simple diagonal matrix with a linear progression from one to 
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forty-two was used. Figure 3.4 shows the control input to the line conditioner for the 
three highest offending harmonics. As in the MSE case, the values are reached 
asymptotically without any overshoot, thus showing the stability of the algorithm. The 
optimal estimation algorithm demonstrated superior robustness and stability compared to 
that of the MSE with respect to the linear predictor constant a. The optimal estimator 
provided stable and consistent results for positive values of a up to one. While 
individual harmonics in the MSE case were highly sensitive, and often grew unstable, 


with changes in a. 
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Figure 3.4 Control Input to Conditioner U, 
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Figure 3.5 shows the total harmonic distortion for several values of a as a function 


of time. 


Total Harmonic Distortion 


0.2 0.4 : 
Time (sec) 


Figure 3.5 Total Harmonic Distortion (THD) 
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An additional piece of information provided by the optimal estimation algorithm 
is an adaptive estimate of the load impedance in the form of a matrix, H. Table 3.1 
gives a break down of several of the actual and estimated matrix values along with their 
respective errors. 


TABLE 3.1 ACTUAL AND ESTIMATED IMPEDANCE MATRIX COEFFICIENTS 


| Actual | 6 {| 7 | 10 | | ta | is | te {9 | 2 | 23 
| Estimated | 6 | 7 | to | | ia [as {is | 19 | 2 | 23 












It is important to point out that only harmonic frequencies which are present in the 
offending noise will produce impedance coefficient estimates for the H matrix. This is 
simply the result of having no error to drive the RLS equation. If noise other than an 
odd harmonic of the fundamental were present an estimate of that impedance would 


appear in the matrix H. 
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IV. CONCLUSIONS 


A. MSE ALGORITHM VS. OPTIMAL ESTIMATOR 

When the load is linear, the optimal estimation algurithm proved to be much more 
effective in eliminating the associated harmonic noise and more robust with respect to 
changes in the gain of the linear prediction, as indicated in Equations 2.8 and 3.10. In 
addition to elimination of the harmonic noise, the optimal estimator provides a linear 
representation of the present system load in the form the impedance matrix H. The 
impedance information of the H matrix is beneficial in helping to determine the proper 


specifications of the active power line conditioner for the particular application. 


B. FUTURE RESEARCH 

Because the optimal estimation algorithm uses the Fast Fourier Transform of the 
error signal it is currently limited to the linear case. Nonlinearities in the error signal 
present a difficult obstacle to overcome using standard Fouier transforms. Investigation 


into adapting the optimal estimation algorithm for nonlinear contingencies would be 


highly beneficial and would provide a better and more accurate model of the power line 


load impedance. 





APPENDIX A - (MATLAB) ADAPTIVE LINEAR MODEL 


THESIS PROGRAM 1 

JOEL ZUPFER 

17 MAY 93 

REVISED i8 NOVEMBER 93 
SIMULATION OF CIRCUIT FIGURE 2.3 


= 700; % MODELLED RESISTOR VALUE (OHMS) 
= 1.857; %MODELLED INDUCTOR VALUE (HENERY’S) 
= 60; % FUNDAMENTAL FREQUENCY (HZ) 
= 120; % NUMBER OF SAMPLES PER CYCLE 
— 120; %LAG FOR CALCULATION OF THE MSE 
= input(’Number of cycles = °); %NUMBER OF CYCLES IN SIMULATION 


input(’Number of cycles between gradient calcu!ations = °); 
%CYCLES BETWEEN GRADIENT CALCULATIONS 


= N*round(Cycle/2)*check; %TOTAL NUMBER OF POINTS IN SIMULATION 
= N*check; % COUNTER FOR CHECKING ERROR VOLTAGE 
= 3e-7; % ACCELERATION FACTOR 


[.51 .16 .056 .035 .025 .025 .02 .02 .015 .015]’; 

®INITIAL WEIGHTS OF ODD HARMONICS 
(3579 {1 13 15 17 19 21}; ODD HARMONICS OF FUNDAMENTAL FREQUENCY 
pi*[1/3 1/4 1/5 3/2 4/2 5/2 1/6 2/6 3/6 8/6)’; %HINITIAL PHASE OF ODD HARMONICS 


a 


H.*(-1). *accel; % ACCELERATION FACTOR ADJUSTED FOR HARMONICS 

% % % % % % Fo % % % % % % KINITIALIZE PARAMETER MEMORY % %%%& %& % % % HHH % % % % % % 
= zeros(length(Wt),totN); % HARMONIC CURRENT 
= zeros(length(Wt),totN); % CONDITIONER CURRENT 
= zeros(1,totN); % FUNDAMENTAL CURRENT 
= zeros(iength(Wt),totN-1); %ERROR VOLTAGE 
= zeros(1,totN-1); %FUNDAMENTAL VOLTAGE 
= zeros(length(Wt),round(Cycle/2) +2); %HARN.UNIC SINE WEIGHTS 
= zeros(length(Wt),round(Cycle/2) + 2); %HARMONIC COSINE WEIGHTS 
= zeros(length(Wt),Cycle+ 1); %MEAN SQUARE ERROR 
= zeros(1,Cycle+ 1); %TOTAL HARMONIC DISTORTION 
= zeros(length(Wt),round(Cycle/2)); %GRADIENT OF SINE HARMONICS 
= zeros(length(Wt),round(Cycle.2)); %GRADIENT OF COSINE HARMONICS 


% % % Yo To Fo Fo % % Fo So K Ko % % % % % KSINITLALCONDITIONS % %& % % % % H% FF % Ho % Io Fo Go Fo Fo Re Fo % 


Sh(:,2) 
Ch(:,2) 


= ones(length(Wt),1)./H.*.1; %INITIAL WEIGHTING FACTOR FOR SINE HARMONIC 
= ones(length(Wt), 1)./H.*. 1%INITIAL WEIGHTING FACTOR FOR COSINE HARMONIC 


If(1) = 10*sqrt(2); 
1f(2) = 9.9*sqrt(2); %PEAK CURRENT VALUES 


31 


% Fo To Fo To To Fo Fo To Fo To To To To %o So To %e MAIN PROGRAM % % % % % % Jo Yo To To Fo To Fo To Fo To Go %o So Mo Fo % 

To %o Fo To To %o %o So Go To Fo BUS VOLTAGE WITH NO CONDITIONER CURRENT % % % % % % & % %% 

fork = 3:N-1 
sample 
Tf(k) 


H.*2*pi*k/N; % DISCRETE SAMPLE POINT 
10*sqrt(2)*cos(2*pi*k/N); % FUNDAMENTAL CURRENT 
Ih(:,k) Wt.*cos(sample + P); %HARMONIC CURRENT OF THE LOAD 
Ic(: ,k) (Sh(:,1).*sin(sample) + Ch(:,1).*cos(sample)); % CONDITIONER CURRENT 
Vf(k-1) = (If(k) - If(k-2))*(L*N*f/2) + If(k-1)*R; %FUNDAMENTAL VOLTAGE 
Ve(:,k-1) = (th(:,k) + Ic(:,k) - Ih(:,k-2) - Ic(:,k-2)).*(L*N*f/2) + (Ih(:,k-1) + Ic(:,k-1)).*R; 

%BUS ERROR VOLTAGE 


| | 


end 


MSE(:,1) sqrt(mean(Ve(:, 1:k-1)’.*2))’; %MEAN SQUARE OF BUS ERROR VOLTAGE 
thd(1) sqrt(sum(max(V(: ,2:k-1)").*2/2)*100/(max(V f(2:k-1))/sqrt(2)); 

%TOTAL HARMONIC DISTORTION 
final= k+1; %STEP INDEX 


% % % So %o Go % So % % % % %BUS VOLTAGE WITH CONDITIONER CURRENT % % % % % % % % % %o % % 
for index = I|:round(Cycle/2) 


for k = final:Delay + final-1 
sample = H.*2*pi*k/N; % DISCRETE SAMPLE POINT 
If(k) 10*sqrt(2)*cos(2*pi*k/N); % FUNDAMENTAL CURRENT 
Th(:,k) Wt.*cos(sample + P); % HARMONIC CURRENT OF THE LOAD 
Ic(:,k) (Sh(: index +1).*sin(sample) + Ch(:,index).*cos(sample)); 
%CONDITIONER CURRENT 
Vf(k-1) (If(k) - If(k-2))*(L*N*f/2) + T€(k-1)*R; %FUNDAMENTAL VOLTAGE 
Ve(.,k-1) = (th(:,k) + Ic(:,k) - Ih(:,k-2) - Ic(:,k-2)).*(L*N*f/2) + (Ih(:,k-1) + Ic(:,k-1)).*R; 
%BUS ERROR VOLTAGE 
end 


MSE(:,2*index) = sqrt(mean(Ve(:,k-1-lag:k-1)’.*2))’; 2 MEAN SQUARE OF ERROR VOLTAGE 
thd(2*index) = sqrt(sum(max(Ve(:,k-1-lag:k-1)’).*2/2)*100/(max(Vf(k-1-lag:k-1))/sqrt(2)); 
%TOTAL HARMONIC DISTORTION 
final = k+1; %STEP INDEX 
GradSh(:,index) = (MSE(:,2*index) - MSE(:,2*index-1))./(Sh(:,index+1) - Sh(:,index)); 
%SIN HARMONICS GRADIENT CALCULATION 
Sh(:,index+2) = Sh(:,index+1) - GradSh(:,index).*mu.*MSE(: ,2*index); 
%PREDICTED SINE WEIGHTING FACTOR 


for k = final:Delay + final-1 
sample H.*2*pi*k/N; % DISCRETE SAMPLE POINT 
If(k) 10*sqrt(2)*cos(2*pi*k/N); %FUNDAMENTAL CURRENT 
Th(: ,k) Wt.*cos(sample + P); %HARMONIC CURRENT OF THE LOAD 
Ic(: .k) (Sh(:,index+1).*sin(sample) + Ch(:,index + 1).*cos(sample)); 
% CONDITIONER CURRENT 
(If(k) - If(k-2))*(L*N*f/2) + Hf(k-1)*R; %FUNDAMENTAL VOLTAGE 
(th(:,k) + Ic(:,k) - Ih(:,k-2) - Ic(:,k-2)).*(L*N*f/2) + (Ih(:,k-1) + Ie(:,k-1)).*R; 
%BUS ERROR VOLTAGE 


Vf(k-1) 
Ve(:,k-1) 


end 

















MSE(:,2*index + 1)= sqrt(mean(Ve(:,k-1-lag:k-1)’.*2))’"%MEAN SQUARE OF ERROR VOLTAGE 
thd(2*index+1) | = sqrt(sum(max(Ve(:,k-1-lag:k-1)’).*2/2)*100/(max(V f(k-1-lag:k-1))/sqrt(2)); 
%TOTAL HARMONIC DISTORTION 
final = k+l; %STEP INDEX 
GradCh(:,index) = (MSE(:,2*index+1) - MSE(:,2*index))./(Ch(:,index+1) - Ch(:,index)); 
% COSINE HARMONICS GRADIENT CALCULATION 
Ch(., index +2) = Ch(:,index+1) - GradCh(: ,index).*mu.*MSE(: ,2*index + 1); 
% PREDICTED COSINE WEIGHTING FACTOR 
end 


% %o To To To To To To To To Yo To Fo To To Yo To To So YoPLOT SECTION % % % % % % % %o %o Jo To %o %o To Yo To To To %o Yo %o % 
plot(Ve’);title( Error Voltage’); 

xlabel(’Samples (N)’);ylabel(’ Voltage (V)’);pause 

plot(Ic’);title(’ Conditioner Current’); 

xlabel(’Samples (N)’);ylabel(’Current (I)’);pause 

plot(MSE’);grid;title(’ Expectation’); grid; 

xlabel(’Samples (N)’);ylabel(’ Magnitude’);pause 

plot(Sh’);title(’Sin weighting coefficients’); grid; 

xlabel(’Samples (N)’);ylabel(’ Magnitude’);pause 

plot(Ch’);title( Cosine weighting coefficients’); grid; 

xlabel(’Samples (N)’);ylabel(’ Magnitude’); 
plot(GradSh’);grid;title(’Grad Sh’); 
xlabel(’Samples (N)’);ylabel(’ Magnitude’ );pause 
plot(GradCh’);grid;title(’Grad Ch’); 
xlabel(’Samples (N)’);ylabel(’ Magnitude’); 
plot(thd);title(’ Total Harmonic Distortion’);grid 
xlabel(’ Number of Cycles’);ylabel(’ Percent (%)’); 


APPENDIX B - (MATLAB) ADAPTIVE NONLINEAR MODEL 


THESIS PROGRAM 2 (NON LINEAR LOAD) 
JOEL ZUPFER 

28 MAY 93 

REVISED 23 NOVEMBER 93 

SIMULATION OF CIRCUIT FIGURE 2.4 


700; %MODELLED RESISTOR VALUE (OHMS) 
1.857; %MODELLED INDUCTOR VALUE (HENERY’S) INITIAL (NONLINEAR L) 
60; % FUNDAMENTAL FREQUENCY (HZ) 
120; % NUMBER OF SAMPLES PER CYCLE 
120; %LAG FOR CALCULATION OF THE MSE 
input(’Number of cycles = °); % NUMBER OF CYCLES IN SIMULATION 
input(’ Number of cycles between gradient calculations = °); 

%CYCLES BETWEEN GRADIENT CALCULATIONS 
N*Cycle*check %TOTAL NUMBER OF POINTS IN SIMULATION 
N*check; %COUNTER FOR CHECKING ERROR VOLTAGE 
= 3e-7; % ACCELERATION FACTOR 
= [.51 .16 .056 .035 .025 .025 .02 .02 .015 .015]’; 

%INITIAL WEIGHTS OF ODD HARMONICS 
= [3579 11 13 15 17 19 21}’; S%ODD HARMONICS OF FUNDAMENTAL FREQUENCY 
= pi*(1/3 1/4 1/5 3/2 4/2 5/2 1/6 2/6 3/6 8/6]’; %INITIAL PHASE OF ODD HARMONICS 
= H.“(-1).*accel; % ACCELERATION FACTOR ADJUSTED FOR HARMONICS 


%o To To To To %o To Yo To To To Fo %o Yo % TAINITIALIZE PARAMETER MEMORY % % % % % % % % % % % % % %o % 
th = zeros(length(Wt),totN); %HARMONIC CURRENT 
zeros(length(Wt),totN); %CONDITIONER CURRENT 
zeros(1,totN); % FUNDAMENTAL CURRENT 
zeros(length(Wt),totN-1); %BUS VOLTAGE 
zeros(length(Wt),round(Cycle/2) +2); %HARMONIC SIN WEIGHTS 
zeros(length(Wt) round(Cycle/2) +2); %HARMONIC COS WEIGHTS 
zeros(length(Wt),Cycle+ 1); %MEAN SQUARE ERROR 
zeros(1,Cycle+ 1); %TOTAL HARMONIC DISTORTION 
zeros(length(Wt),Cycle); %GRADIENT OF SIN HARMONICS 
zeros(length(Wt),Cycle); %GRADIENT OF COS HARMONICS 


iouu wv on 


%o To %o To Yo So %o To %o %o %o %o To To % % To %o TINITIALCONDITIONS% % % % % % % % % %% % % % % % % % %% 
Sh(:,2) = ones(length(Wt),1)./H*.1; | %INITIAL WEIGHTING FACTOR FOR SINE HARMONICS 
Ch(:,2) = ones(length(Wt), 1)./H*.1-%INITIAL WEIGHTING FACTOR FOR COSINE HARMONICS 
If(1) 10*sqrt(2); 

If(2) = 9.9*sqrt(2); PEAK CURRENT VALUES 











To To Fo To To Fo Fo To %o To Fo Fo Fo To To Fo To % FMAIN PROGRAM % % % % % To To Fo Fo Fo Fo Fo Fo Fo To Yo Fo Fo Fo Ro % % 


% To %o To %o %o So %o % So % %BUS VOLTAGE WITH NO CONDITIONER CURRENTS % % %%% %%HEM 
fork = 3:N-1 


sample = H.*2*pi*k/N; % DISCRETE SAMPLE POINT 
If(k) = 10*sqrt(2)*cos(2*ni*k/N); % FUNDAMENTAL CURRENT 
Th(:,k) = Wt.*cos(sample + P) %HARMONIC CURRENT OF THE LOAD 
Ic(:,k) = Sh(:,1).*sin(sample) + Ch(:,1).*cos(sample); % CONDITIONER CURRENT 
L = Li *ones(10,1) ./ (1 +abs((Ih(:,k-1) +Ic(:,k-1))/10)); % NONLINEAR INDUCTANCE 
Vf(k-1) 9 = (If(k) - If(k-2))*(Li*N*£/2) + Hf(k-1)*R; % FUNDAMENTAL VOLTAGE 


Ve(:,k-1) = (th(:,k) + Ice(:,k) - Ih(:,k-2) - Ic(:,k-2)).*(L*N*f/2) + (ih(:,k-1) + Ic(:,k-1)).*R; 
%BUS VOLTAGE 
end 


MSE(:,1)= sqrt(mean(Ve(:,1:k-1)’.*2))’; %MEAN SQUARE OF ERROR VOLTAGE 
thd(1) = sqrt(sum(max(Ve(:,2:k-1)’).*2/2)*100/(max(Vf(2:k-1))/sqrt(2)); 

%TOTAL HARMONIC DISTORTION 
final = k+l; %STEP INDEX 


% To Go %o To To To To To To So %o KBUS VOLTAGE WITH CONDITIONER CURRENT % % % % % % % % % %o % % 
for index = 2:Cycle 


for k = final:Delay +final-1 


sample = H.*2*pi*k/N; % DISCRETE SAMPLE POINT 
If(k) = 10*sqrt(2)*cos(2*pi*k/N); %FUNDAMENTAL CURRENT 
Th(:,k) = Wt.*cos(sample + P); %HARMONIC CURRENT OF THE LOAD 
Ic(: ,k) = Sh(:,index).*sin(sample) + Ch(:,index-1).*cos(sample); 

% CONDITIONER CURRENT 
L =Li *ones(10,1) ./ (1 +abs((Ih(:,k-1) +Ic(:,k-1))/10)% NONLINEAR INDUCTANCE 
Vftk-1) = (If(k) - If(k-2))*(Li*N*f/2) + If(k-1)*R; %FUNDAMENTAL VOLTAGE 


Ve(:,k-1) = (Ih(:,k) + Ic(:,k) - Ih(:,k-2) - Ic(:,k-2)).*(L*N*f/2) + (Ih(:,k-1) + Ic(:,k-1)).*R; 
%BUS VOLTAGE 
end 
MSE(:,2*index-2)= sqrt(mean(Ve(: ,k-1-I:k-1)’.*2))’; 
%MEAN SQUARE OF ERROR VOLTAGE 


thd(2*index) = sqrt(sum(max(Ve(: ,k-1-lag:k-1)’).*2/2)*100/(max(V f(k-1-lag:k-1))/sqrt(2)); 
%TOTAL HARMONIC DISTORTION 
final = k+l; %STEP INDEX 


GradSh(:,index) = (MSE(:,2*index-2) - MSE(:,2*index-3))./(Sh(:,index) - Sh(:,index-1)); 
%SINE GRADIENT CALCULATION 
Sh(:,index+1) | = Sh(:,index) - GradSh(:,index).*mu.*MSE(:,2*index-2); 
%PREDICTED SINE WEIGHTING FACTOR 
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for k = final:Delay + final-1 


sample = H.*2*pi*k/N; % DISCRETE SAMPLE POINT 
If(k) = 10*sqrt(2)*cos(2*pi*k/N); %FUNDAMENTAL CURRENT 
Ih: ,k) = Wt.*cos(sample + P); %HARMONIC CURRENT OF THE LOAD 
Ie(:,k) = Sh(:,index).*sin(sample) + Ch(-,index).*cos(sample); 
% CONDITIONER CURRENT 
L = Li *ones(10,1) ./ (1 +abs((Ih(: ,k-1)+Ic(:,k-1))/10)## NONLINEAR INDUCTANCE 
Vf(k-1) = = (If(k) - If(k-2))*(Li*N*f/2) + If(k-1)*R; %FUNDAMENTAL VOLTAGE 
Ve(:,k-1) = (fh(:,k) + Ic(:,k) - Ih(:,k-2) - Ic(:,k-2)).*(L*N*f/2) + (Ih(:,k-1) + Ic(:,k-1)).*R; 
%BUS VOLTAGE 
end 


MSE(:,2*index-1)= sqrt(mean(Ve(: ,k-1-1:k-1)’.*2))’; 
%MEAN SQUARE OF ERROR VOLTAGE 
thd(2*index + 1)= sqrt(sum(max(Ve(:,k-1-lag:k-1)’).*2/2)*100/(max(Vf(k- 1 -lag:k-1))/sqrt(2)); 
%TOTAL HARMONIC DISTORTION 
final = k+1; %STEP INDEX 
GradCh(:,index) = (MSE(:,2*index-1) - MSE(:,2*index-2))./(Ch(:,index) - Ch(:,index-1)); 
% COSINE GRADIENT CALCULATION 
Ch(:,index+1)= Ch(:,index) - GradCh(:,index).*mu.*MSE(:,2*index-1); 
%PREDICTED COSINE WEIGHTING FACTOR 
end 


% % % To %o To To To To To To So To To To So To % YPLOTTING SECTION % % % % % %o % %o %o %o To To To %o To Yo To To Fo % 
plot(Ve’);title(’ Error Voltage’); 

xlabel(’Samples (N)’);ylabel(’ Voltage (V)’);pause 
plot(Ic’);title(’ Conditioner Current’); 
xlabel(’Samples (N)’);ylabel(’Current (I)’);pause 
plot(MSE’);grid;title(’ Mean Square Error’); grid; 
xlabel(’Samples (N)’);ylabel(’ Magnitude’);pause 
plot(Sh’);title(’Sin weighting coefficients’); grid; 
xlabel(’Samples (N)’);ylabel(’ Magnitude’);pause 
plot(Ch’);title(’Cosine weighting coefficients’); grid; 
xlabel(’Samples (N)’);ylabel(’ Magnitude’); 
plot(GradSh’); grid;title(’Grad Sh’); 

xlabel(’Samples (N)’);ylabel(’ Magnitude’);pause 
plot(GradCh’);grid;title(’Grad Ch’); 
xlabel(’Samples (N)’);ylabel(’ Magnitude’); 
plot(thd);title( Total Harmonic Distortion’); grid 
xlabel(’Number of Cycles’);ylabel(’Percent (%)’); 
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APPENDIX C - (MATLAB) OPTIMAL ESTIMATION MODEL 


%% THESIS PROGRAM 3 

%% JOEL ZUPFER 

%% 1 JUNE 93 

%% REVISED 30 NOVEMBER 93 

%% SIMULATION OF CIRCUIT FIGURE 2.3 


clear 
R = 700; %MODELED RESISTOR VALUE (OHMS) 
L = 1,857; %MODELED INDUCTOR VALUE (HENERY'’S) 
f = 60; %FUNDAMENTAL FREQUENCY (HZ) 
wt = (00 .51 0 .16 0 .056 0 .035 0 .025 0 .025 0 .02 0 .02 0 .015 0 .015]; 
% WEIGHTS OF THE HARMONICS INITIAL CONDITION 
H = {123456789 1011 12 13 14 15 16 17 18 19 20 21); 
%HARMONICS OF FUNDAMENTAL FREQUENCY 
P = pi*(0 0 1/3 0 1/40 1/5 0 3/2 0 4/2 0 5/2 0 1/6 0 2/6 0 3/6 O 8/6); 
% PHASE OF HARMONICS INITIAL CONDITION 
M = 2*length(H); %NUMBER OF SAMPLE POINTS/PERIOD 
ALPHA = 0.75; % WEIGHTING FACTOR FOR CONTROL DETERMINATION 
Cycle = input(’Number of cycles = '); %NUMBER OF CYCLES IN SIMULATION 


% To % %o To To To To To To To To Yo Mo TANITIALIZE PARAMETER MEMORY % % % % % % %o % % % %o % % % 


D = [1:42]; % DIAGONAL ELEMENTS OF H 
HL = diag(D); % TRANSFER MATRIX OF SYSTEM LOAD 
Qhat = eye(M)*.1; %INITIAL GUESS FOR INVERSE OF H 
Q = zeros(2*length(H),2*length(H)); %INITIAL VALUE OF Q MATRIX 
theta0 = zeros(1,1); %THE ESTIMATION FOR THE DC COMPONENT 
theta42 = zeros(1,1); %THE ESTIMATION FOR THE 21st HARMONIC 
theta = zeros(2,2*(length(H)-1)); %THE ESTIMATION OF THE SYSTEM IMPEDANCE 
GO = 50; %INITIAL VALUE OF GAIN MATRIX(MEASURE OF UNCERTAINTY) 
G42 = 50; %INITIAL VALUE OF GAIN MATRIX 
G = [ones(2*H(length(H)-1), 1),zeros(2*H(length(H)-1),2),ones(2*H (length(H)-1), 1)]*50; 
%INITIAL VALUE OF GAIN MATRIX 
If = zeros(1,M); %INITIALIZE THE FUNDAMENTAL CURRENT 
Vf = zeros(1,M); %INITIALIZE THE FUNDAMENTAL VOLTAGE 
Ve = zeros(1,M); %INITIALIZE THE NOISE VECTOR 
xX = zeros(M,M); %INITIALIZE THE MATRIX OF SUM OF COS & SIN 
E = zeros(M,Cycle); % ACTUAL ERROR 
Ek = zeros(M,Cycle); % PERIODIC ERROR 
Ehat = zeros(M,Cycie); % ESTIMATED PERIODIC ERROR 
Uk = zeros(M,Cycle); %PERIODIC CONTROL INPUT 
Uhat = zeros(M,Cycle); % ESTIMATED PERIODIC CONTROL INPUT 
Yk = zeros(M,Cycle); %PERIODIC SYSTEM OUTPUT 
thd = zeros(1,Cycle); %TOTAL HARMONIC DISTORTION 
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HKD VRHHEHHHKHNGKSHKRMAIN PROGRAMRKEEEERERKEEKREEEEES 
for n = 0:M-1 


If(n+1) = 10*sqrt(2)*cos(2*pi*n); 
Sample = H.*(2*pi*n/M); 
X(:,n+1) = reshape({cos(Sample);sin(Sample)},2*length(H), 1); 
%MATRIX FOR SAMPLED VALUE AS SUM OF COS & SIN 
end 
WwW = reshape([cos(P).*Wt;sin(P). *Wt],2*length(H), 1); 

%HARMONIC CURRENT WEIGHT VECTOR 
th = W’'*X; %THE HARMONIC CURRENT OF ONE PERIOD 
Vf(l) = = (If(2) - If(M))*(L*M*f/2) + If(1)*R; %FUNDAMENTAL VOLTAGE 
Ve(:,1) = (th(:,2) - Th(:,M)).*(L*M*f/2) + th(:,1).*R; %ERROR VOLTAGE 
forn = 3:M 

Vf(n-1) == (If(n) -If(n-2))*(L*M*f/2) + If(n-1)*R; %FUNDAMENTAL VOLTAGE 
Ve(:,n-1) = (fh(:,n) - Ih(:,n-2)).*(L*M*f/2) + th(:,n-1).*R; %ERROR VOLTAGE 
end 
Vf(M) = = (If(1) - If(M-1))*(L*M"*f/2) + If(M)*R; %FUNDAMENTAL VOLTAGE 
Ve(:,M) = (ih(:,1) - Th(:,M-1)).*(L*M*f/2) + Th(:,M).*R; %ERROR VOLTAGE 
MAGf = max(Vf)/SQRT(2); %RMS VALUE OF THE FUNDAMENTAL 
Yk(:,1) = H*UKkK(:,1); %SYSTEM OUTPUT 
E(:,1) = Ve’; % ACTUAL ERROR 
Ek(:,1) = (fftrig(E(:,1)))’; %PERIODIC ERROR COEFFICIENTS FREQUENCY DOMAIN 
Uk(:,2) = Uk(:,1) - ALPHA*Qhat*Ek(:,1); % CONTROL COEFFICIENT UPDATE 
for I = L:length(Ek(:,1))/2 -1 
MAG(i,1) = sqrt(sum(Ek(2*i:2*i + 1,1)).*2)); %SUM OF THE ERROR VOLTAGES 
end 
MAG(length(Ek(:,1))/2,1) = Ek(ength(Ek(:,1)),1); 
thd(1) = sqrt(sum(MAG.“2/2))*100/(MAGf); %TOTAL HARMONIC DISTORTION 
for k = 2:Cycle %RECURSIVE LOOP FOR ERROR CALCULATION 
Ek(:,k) = fftrig(Ve)’ + H*Uk(:,k); % ACTUAL ERROR IN FREQUENCY DOMAIN 
E(k) = = ifftrig(Ek(:,k))’; %TIME DOMAIN OF ERROR 
Ehat(:,k) = Ek(:,k) - Ek(:,k-1); % ESTIMATED ERROR 
Uhat(:,k)= Uk(:,k) - Uk(:,k-1); %ESTIMATED CONTROL 


thetaQ = theta0 + (Uhat(1,k) - Ehat(1,k)*theta0)*GO*Ehat(1,k)/(1 + Ehat(1,k)*GO*Ehat(1,k)); 
%RECURSIVE LEAST SQUARES ESTIMATION OF COEFFICIENT 

GO = GO - G0*Ehat(1,k)*Ehat(1,k)*GO/(1 + Ehat(1,k)*GO*Ehat(1,k)); 
%UPDATE OF THE GAIN 


Qdi,1) = thetad; 
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end 


forh = 1:2*H(length(H)-1) 


gain = reshape(G(h,:),2,2); %GAIN MATRIX IS 2X2 
if rem(h,2) ==1; %FOR THE ODD LOOPS 
Phi = Ehat(h+ !:h+2,k)’.*[1 -1]; 
else 
Phi = flipud(Ehat(h:h+ 1,k))’; %FOR THE EVEN LOOPS 
end 


theta(:,h) = theta(:,h) + (Uhat(h+1,k) - Phi*theta(:,h))*gain*Phi’/(1 + Phi*gain*Phi’); 
% RECURSIVE LEAST SQUARES ESTIMATION OF COEFFICIENTS 


gain = gain - gain*Phi’*Phi*gain/(1 + Phi*gain*Phi'); % UPDATE OF THE GAIN 
G(h,:) = reshape(gain, | ,4); %CHANGE GAIN MATRIX BACK TO 1X4 
Q(ht+1,h+1:h+2) = theta(:,h)’; % PLACE ESTIMATES WITH IN THE Q MATRIX 
end 
theta42 =theta42 + (Uhat(42,k) - Ehat(42,k)*theta42)*GO*Ehat(42,k)/... 
(1+ Ehat(42,k)*G0O*Ehat(42,k)); 
% RECURSIVE LEAST SQUARES ESTIMATION OF COEFFICIENT 
G42 = G42 - G42*Ehat(42,k)*Ehat(42,k)*G42/(1 + Ehat(42,k)*G42*Ehat(42,k)); 
% UPDATE OF THE GAIN 
Q(42,42) = theta42; % PLACE ESTIAMTES WITH IN THE Q MATRIX 
Uk(:,k+1) = Uk(:,k) - ALPHA*Q*EK(:,k); %UPDATE THE CONTROL 
Yk(:,k) = ifftrig(H*Uk(:,k))’; %SYSTEM OUTPUT 
Nfk = fftrig(Ve + Yk(:,k)’)’; %ERROR AFTER NEW CONTROL INPUT 


fori = 1:length(Nfk(:,k))/2 - 1 


end 


MAG(i,k) = sqrt(sum(Nfk(2*i:2*i + 1),k).*2)); % MAGNITUDE OF ERROR VECTOR 


MAG (length(Nfk(:,k))/2,k) = Nfk(length(Nfk(:,k)),k); %SUM OF THE ERROR VOLTAGES 
thd(k) = sqrt(sum(MAG(-,k).*2/2))*100/(MAGf); %TOTAL HARMONIC DISTORTION 


%o To To To To To To Fo Fo %o Fo % %o Fo Go Go Fo %oe VPLOT SECTIONRRRKLKRKHKKNNMDND LG % GK % % % % 
plot(E’);title(’ Actual Error’);grid; 

xlabel(’Samples (N)’);ylabel(’ Voltage (V)’);pause 
plot(Ek’);title(’ Periodic Error’);grid; 

xlabel(’Samples (N)’);ylabel(’ Voltage (V)’);pause 
plot(Ehat’);title(’ Estimated Error’);grid 

xlabel(’Samples (N)’);ylabel(’ Voltage (V)’);pause 
plot(Uk’);title(’Control weighting coefficients’); grid; 
xlabel(’Samples (N)’);ylabel(’ Magnitude’);pause 
plot(Uhat’);title( Estimated Control weighting coefficients’); grid; 
xlabel(’Samples (N)’);ylabel(’ Magnitude’);pause 
plot(Yk’);title(’System output’); grid; 

xlabel(’Samples (N)’);ylabel(’ Magnitude’); 

plot(thd); title(’ Total Harmonic Distortion’);grid 

xlabel(’ Number of Periods (k)’);ylabel(’ Percent (%)’); 
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(function w = fftrig(x) 

% w=fftrig(x) 

% it computes the fft coefficients of the vector x 
% in trigonometric form. 

% x(n)=wi +w2*cos(2pi/N) + w3*sin(2pi/N) + v4*cos(4pi/N) + w5*sin(4pi/N) +... 
% ... + w(N-2)cos(2pi(N/2-1)/N)+w(N-1)sin(2pi(N/2-1)/N)+w(N)cos(pi 1); 

% where N=length(x), assumed to be a power of 2 (if not x is padded with 0’s) 


X  =fft(x); 

N = =length(X); 
Xr  =real(X); 
Xi =imag(X); 


w(l) =Xr(1)/N; 

k=1:1:(N/2)-1; 
w(2*k) =2*Xr(2:(N/2))/N; 
w(2*k+1) =-2*Xi(2:(N/2))/N; 
w(N) =Xr((N/2)+1)/N; 

end % fftrig 


function x=ifftrig(w) 

% x=ifftrig(w) 

% compute time sequence x from trig. coefficients w. 
% See FFTRIG 


N =length(w); 
Xr(1) =w(1)*N; 
Xi(l) =0; 


k=1:1:(N/2)-1; 
Xr(2:(N/2)) =N*w(2*k)/2; 
Xi(2:(N/2))  =-N*w(2*k+1)/2; 
Xr((N/2)+1) =N*w(N); 
Xi((N/2)+1) =0; 


XxX =Xr+sqrt(-1)*Xi; 

X(N-k+1) = =conj(X(k+1)); 

x =real(ifft(X)); 
end 4% ifftrig 
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